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Random Variables 

Ken Umeno, Aki-Hiro Sato 



Abstract 



This study proposes a pseudo random number generator of (/-Gaussian random variables for a range 



of q values, — oo < q < 3, based on deterministic chaotic map dynamics. Our method consists of 
chaotic maps on the unit circle and map dynamics based on the piecewise linear map. We perform the 
(/-Gaussian random number generator for several values of q and conduct both Kolmogorov-Smirnov 



■ (KS) and Anderson-Darling (AD) tests. The (/-Gaussian samples generated by our proposed method pass 

Q 

the KS test at more than 5% significance level for values of q ranging from -1.0 to 2.7, while they pass 

' the AD test at more than 5% significance level for q ranging from -1 to 2.4. 
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Index Terms 



ly-j ' Map dynamics, Chebyshev polynomials, pseudo random number generator, (/-Gaussian distribution, 

| ergodic theory. 
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j^-j , I. Introduction 

rS 

The (/-Gaussian distributions have been studied in a wide variety of fields from natural sciences to social 
sciences. They have been applied in thermodynamics, biology, economics, and quantum mechanics. The 
generating mechanism is still an open question, but several mechanisms that have been shown to produce 
(/-Gaussian distributions are known, such as multiplicative noise, weakly chaotic dynamics, correlated 
anomalous diffusion, preferential growth of networks, and asymptotically scale-invariant correlations flT|. 
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In the heavy-tail domain (1 < q < 3), the ^-Gaussian distribution is equivalent to the Student's t- 
distribution. 

In the context of finance, the (/-Gaussian distribution (1 < g < 3) is referred to as a Student's t 
distribution J2). This is commonly used in finance and risk management, particularly to model conditional 
asset returns of which the tails are wider than those of normal distribution. The distribution is also known 
as Pearson Type-II (for compact support (q < 1) and Type VII (infinite support (q > 1) |3 |. For example, 
Bollerslev used the Student's t to model the distribution of the foreign exchange rate returns 0. Bening 
and Korolev provide an instance where the distribution is appropriate as a model, i.e. the case of random 
sample sizes (51 . Vignat and Plastino obtained similar results 0. Other work attempts to show the q- 
Gaussian distribution as an attractor in the context of dependent systems [7]. Moreover, Umarov et al. 
consider a (/-extension of ct-stable Levy distribution JSJ. 

More recently, (/-Gaussian distributions have been derived from the maximization of non-extensive 
entropy 0]] and studied in the context of the generalization of Gauss' Law of Errors ||9). (/-Gaussian 
distributions can be derived from an infinite normal mixture with an inverse gamma distribution. This 
concept is known as superstatistics in non-equilibrium thermodynamics ifTOI . (/-Gaussian distributions 
also appear as unconditional distributions of multiplicative stochastic differential equations |[TTI . 

Recently, the generalized Box-Muller method (GBMM) was proposed by Thisleton et al. lPT3l . Their 
method uses transformation including the (/-logarithmic, sine, and cosine functions in terms of uniform 
random variables. Here, based on the ergodic theory lfl4l of dynamical systems, we propose a family of 
chaotic maps with an ergodic invariant measure given by (/-Gaussian density. Ulam and von Neumann 
considered the logistic map X n+ \ = 4X n (l — X n ) in the late 1940s, and found its randomness ITT31 . One 
of the authors (K. Umeno) proposed chaotic mechanism to generate power-law random variables iTTol . 
This method can generate power-law random variables in the Levy stable regime from the superposition 
of the random variables. One of the authors (A.-H. Sato) also proposed multiplicative random processes 
to generate power-law random variables iTTTl . Currently, we can use the map dynamics to design random 
sequences with an explicit ergodic invariant measure more precisely lfl8l . lTT9l . 

In this article, we propose a method to generate (/-Gaussian random variables based on deterministic 
map dynamics. Our method is based on ergodic transformations on the unit circle and a map composed of 
the piecewise linear map with both the (/-exponential and (/-logarithmic function. This method is a direct 
method different from Ref. |[T3l . |[T6l and can generate (/-Gaussian random variables for — oo < q < 3 
including infinite variance and infinite mean regimes. We generate (/-Gaussian random variables for several 
cases of q, and conduct statistical testing by means of analytical cumulative distribution functions. 
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II. Review of the generalized Box-Muller Method 
The zero-mean normal (/-Gaussian distribution parameterized by q is described as 

(?<1) 



9(x;q) 



B 



2-g 1 
l-q'2 



^exp-f (q = l) , (1) 



B 

where B(a,b) is the beta function, which is defined as 



1 1 1 
1—q 2'2 



B(a,b) = [ t^n-tf^dt. 
Jo 



(2) 



For q < 1 symmetric distributions with compact support ranging from — J to y' appear. Specif- 
ically, the normalized Wigner distribution is obtained at q = —1. In the case of 1 < q < 3, Equation 
[T]has heavy-tails and g(x;q) « const .\x\ u ~ 1 , where v = (3 — q)/(q — 1) > is related to the degree 
of freedom of the Student's t-distribution. v is coincident with the tail index of the complementary 
cumulative distribution of g(x; q). This also gives an existence condition in the heavy-tail regime of the 
q-Gaussian distribution. 

Firstly, let us start our discussion from the GBMM proposed by Thistleton et al. Ifl3l . To introduce 
their method to generate g-Gaussian random variable, we define a (/-analog of both exponential and 
logarithmic function. 

Definition 1. Suppose the one-dimensional ordinary differential equation 

Ah 

— = h«, h(0) = l. (3) 

The solution is given as 

f (l + (l-g)™)^ l + (l- q ) w>0 
h[w) = \ ■ (4) 

I elsewhere 

We call the solution h(w) (/-exponential function. Obviously, one has 

i 

lim(l + (1 - q)w) ^ = e w . (5) 



Definition 2. We define the inverse function of Equation [4] 

ln» = - " - 1 (w>0), (6) 
1 - q 
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which we call the q-logarifhmic function. Clearly, we get 

. w l ~ q - 1 
lim 



g-+l 1 



lnu>. 



(V) 



Definition 3. The GBMM [13] is given by transformations from i.i.d. uniform random variables u\ 
and U2 ranging from to 1. 



x 

y 



-2 ln q (ui) cos(27TU2) 



-2 ln g (m) sin(2-7rn2) 



Proposition 1. The joint probability density of x and y in Equation [U is given by 

1 



Vxy{x,y) = — exp 2 _i/ g (-|(x 2 +y 2 )). 



(8) 



(9) 



Proof of Proposition 1. From Equation [U we obtain 

Pu x (ux)p U2 (u 2 ) = px,y(x,y) 



d(x,y) 



d(ui,u 2 ) 
2irp x ,Y(x,y)ui Q 



Px,y(x,y) 



1 



1 



exp 9 (--(s +y )) 



2?r 

^exp 2 _ 1/g (--(x +y . 



(10) 



where we used the equality 

{exp q (x)) q = (1 + (1 - = exp 2 _ 1/q (qx). 

Note that Equation [10] is recognized as a two-dimensional (/-normal distribution, 

PxA^v) = h eXPr( C 2 + D\l-r) ix2 + y2) )> 
where r = 2 — l/q and D = 2. This is properly parameterized with each marginal g-variance equal to 



(11) 



(12) 



one. 



Proposition 2. The marginal distribution of x is given by 



Px(x) = < 



1 exp(-^ 



l-^4x 2 

3-q' 



<Vf^ («'<i) 



2tt 



9 

'= 1 



(13) 



1-9' 



(K 4 < 3) 
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where q' = (3q — l)/(g+ 1). Hence, x in Equation [8] gives a (/'-Gaussian random variable. 

Proof of Proposition 2. Integrating Equation [9] in terms of y, we obtain Equation [[3] In the case of 
q = 1, we obviously obtain 



Px{x) = / p x .y{x,y)dy 

J -co 



In the case of 1 < q < 3, we have 

/■oo 

Px(x) = / px,Y{x,y)dy 



oo 



1 Z" 00 / *7 \ 

(c = i+ 2 



?"1 2 



7rV2(g-l) 

Using the equality among beta function and gamma functions 



1 1 I 1 Z" 1 1 _ 3 / ,1 

7rv2(<7 — 1) JO 

0-11 \q-\ 2'2/V 2 / 



r(a + 6) 

where the gamma function is defined as 



/"OO 

r(o) = / t a_1 e-*dt, (17) 
Jo 



and 7T = T(i) 2 , 



7T 



r(i) 2 r(-^) "rdM-aj)- 



(18) 



January 10, 2013 



DRAFT 



6 



Since we have T(^_) = - 1)1^ - 1) = - 1), we get 

q-1 



Therefore, Equation [151 can be rewritten as 



Setting g = we obtain 



-1,|)V 2 M 1+ 2 



,r" 



2\ 2(,-l) 



B 



-1 

i _ I l) \3-q'J 



l-q' 2> 2 



1 + 



1-9' 



(19) 



(20) 



(21) 



In the case of q < 1, we obtain the joint density px,y{x,y) has a compact support ranging from 



2 ^2 



l- 9 



x z to 



2 „2 



1-9" 



Px(x) 



_Px,Y{x,y)dy 



l /V^V l 



7T Jo 



(c 



(l-^^V + I/ 2 )) 1 "'^ 



n 1 - 9 2 
1 x 



1 , 1 f 00 , 3 1 

/Q , n C~ + ^ / *i=i"(l-*)"dt 
7rV2(g-l) Vi 



H) 



1 



1 Dh 1 



2 \ 2(1-,) 



7rV2(g-l) Vl-g'2M- 2 
Similarly to the case of 1 < q < 3 setting g = §£^r, we obtain 

px{x) 



B(u^r- q,) 1 *- q ' 



l-q' 



where \x\ < 



3-q' 
l-q'- 



(22) 



(23) 



Figure \T\ shows the distribution of Equation [9] for several cases of q. The distribution is the spinning 
object. The marginal distribution in terms of £ is also equivalent to Equation [T3l 
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III. Map dynamics 

Adler and Rivlin considered Pd{w) = cos dO, where w = cos 9, < 9 < it, defined by Chebyshev 
polynomial of degree d 11201 , where d is an integer. Clearly, Pd(w) is permutable Pd 1 (Pd 2 (w)) = Pd^ i w )- 
They proved that for d > 2 the ergodic invariant measure of the map dynamics w n+ \ = Pd(w n ) has an 
explicit density function invariant measure Hw{ w ) = T ^/ 1 1 _ tt|2 ■ 

More generally, let us extend the Chebyshev polynomial to a two-dimensional case as lf2TTl . 

Definition 4. Pd{w) and Qd{w,v) are given as real and imaginary parts of binomial expansion, 

(w + iv) d = P d {w)+iQ d {w,v), (24) 

where the equality w 2 + v 2 = 1 is necessary in order to obtain Pd{w) from this expansion. Here, in this 
definition we used the Eular equality 

(exp(i6>)) d = exp(id9) = cos dO + i sin dO. (25) 



This Pd(w) is the Chebyshev polynomial of degree d. The first few polynomials are explicitly given 

by P±(w) = w, Qi(w,v) = v, P2{w) = 2w 2 — 1, Q2(w,v) = 2wv, Ps(w) = 4w 3 — 3w, Qs(w,v) = 
v{4w 2 - 1), P 4 (w) = 8w 4 - 8w 2 + 1, Qi{w,v) = v(8w 3 - Aw), P 5 (w) = 16w 5 - 20w 3 + 5w, 
Q 5 (w, v) = v{\Qw A - 12w 2 + 1), P 6 (w) = 32w 6 - 48w 4 + 18w 2 - 1, Q 6 (w, v) = v(32w 5 - 32w 3 + 6w), 
P 7 (w) = Uw 7 - ll2w 5 + 56w 3 - 7w, Q 7 (w,v) = v(64w 6 - 80w 4 + 24w 2 - 1), P 8 (w) = 128w 8 - 
256w e + 160w A - 32w 2 + 1, and Qs(w, v) = v{\28w 7 - 192w 5 + 80u> 3 - 8w). 
Definition 5. For d > 2, we define the map dynamics 

w n+1 = P d {w n ) (26) 

V n +1 = Qd(Wn,V n ) (27) 

on the unit circle w 2 l + v 2 l = 1. The set of variables (w n ,v n ) is uniformly distributed on the unit circle 
if we set an initial condition of (wq, vq) on the unit circle. We set vq as an arbitrary value in (0, 1) and 



wo is given by = ±y 1 — Vq. 

Lemma 1. The joint density of ergodic invariant measure for u> n +i = Pd{w n ) and v n +i = Qd(w n , v n ) 
follows 

p wy (w,v) = — — =— . (28) 

Ztt\/w z + v A 
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Proof of Lemma 1. In addition to Pd{w) = cos dO, we introduce Qd(w, v) = sin dd, where w = cos 8 
and v = sin 8, < 8 < 2ir. From the equality given in Equation |25l the angular 8 n of w n + iv n follows 
the map dynamics 



>n+ 



i = dd n mod 2ir, 



(29) 



which is ergodic and has an ergodic density function frill 



Pe(d) = — (0<e<2Tr). 
2tt 



(30) 



Transforming the orthogonal coordinate (w,v) into the polar coordinate (a, 8) by w = a cos 8 and 



v = asind, we have Pa{o) = S{a — 1). Since ons has a = \Jw 2 + v 2 ^ = cos 8 



cos!', 



dv 



q- — sin 8, and || = a cos 0, the Jacobian matrix is expressed as 

d{8, a) 



d(w, v) 



costf —asm; 
sin 8 a cos 8 
-l 1 



\A« 2 + 



Therefore, the joint density of the ergodic invariant measure of w and v can be described as 



p W y(w,v) = pe(8)pA(a 



d(8,a) 



d(w, v) 
6(Vw 2 + v 2 - 1) 
2ny/w 2 + f 2 



-a sin f, 



(31) 



(32) 



Lemma 2. The density functions of ergodic invariant measure of Equation |26] and Equation 
respectively, have the form: 



Vw(w) = — u = , (33) 



w 



M") = n 2 - (34) 



Proof of Lemma 2. From Equation |32l we can calculate pw{w) and pv{v) as the marginal distribution 
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in terms of w and v. Integrating Equation [32] with respect to v and w, we respectively obtain 



Pw(w) 



p v {v) 



Pwv(w, v)dv 



-oo 
oo 



S(Vw 2 + v 2 -l) 
2ir\/w 2 + v 2 



dv 



TTa/T 



Pwv{w, v)dw 



-oo 
oo 



5(V 



1) 



2ir\/w 2 + v 2 



dw 



(35) 



(36) 



Definition 6. As an alternative method for generating (/-Gaussian random variables, we propose chaotic 
maps based on the following map dynamics: 



Wn+l 
V n +l 
Zn+l 



Pd(Wn) 

Qd(Wn,V r , 

fl,c{Zn) 



(37) 



where 



fl,c{z) 



wl + vl 



c 

1, z n > 



(38) 
(39) 



assuming 



exp, 



-21n 9 (u), 



(40) 
(41) 
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where 7] (u) is an l-th order piecewise linear map denned as 

Ix (0 < u < l/l) 



Ti{u) 



-lx + 2 



{l/l << 2/1) 



(I : odd) 



Ix-(l-l) (l-l/l <x <1) 
Ix (0 < u < l/l) 

-lx + 2 (l/l << 2/1) 



(42) 



(I : even) 



-lx + l (l-l/l <x <l) 
For example, in the case of / = 2, Equation @2] gives the tent map 

2u (0 < u < 1/2) 

-2n + 2 (1/2<-u<1) 
1- II- 2u\. 



T 2 {u) 



(43) 



In the case of I = 3, Equation 02] is expressed as 

3u (0 < u < 1/2) 

T 3 (it) = < -3n + 2 (1/2 < u < 1) (44) 
3u-2 (2/3<it<l) 

The number of iteration c is an integer greater than or equal to 1. The order / of the piecewise linear 
map is an integer greater than or equal to 2. By using the product among z n , w n , and v n , 

in 



Zjy r i, 



(45) 



we can also obtain two-dimensional deterministic dynamics. The random seed of this pseudo random 
generator is given by (vq, zq), where we set wo as wq = ±^/l — Vq. 

Note that factor 2 in front of ^-exponential function in Equation [43] should be replaced with a value 
both smaller than and close to 2, such as 1.99999, for the round error correction in the case of actual 
numerical computation. 

Lemma 3. The density of ergodic invariant measure of z n+ \ = fi :C (z n ) follows the one-side distribu- 
tion, 

Pz (z) = zexp 2 _ 1/9 (-|z 2 ) (*>0). (46) 
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Proof of Lemma 3. The density of the ergodic invariant measure frPZl of the piecewise linear map 

u n+ i = Ti o • • • o Ti (u n ), (47) 

c 

follows the uniform distribution pu{u) = 1 (0 < u < 1) independently of I and c. Since we obtain 



du/dz = u q ^j — 2 ln q (u) from the transformation in Equation |40] we have 

, , du 



Pz{z) 



PU(Z 



dz 

(exp q (-z 2 /2) 



zexp 2 _ 1/g (--2: 2 ) {z>0). 



(48) 



In this derivation, we used the equality introduced in Equation ITTI 

Theorem 1. The joint density PE^i^f]) °f ergodic invariant measure of map dynamics Equation | 
is the g-Gaussian distribution which is the same as Equation [9] and given by 



Ps,e(£,v) = 7^exp 2 _ 1/(2 (-^(£ 2 + t/ 2 )). 



(49) 



Proof of Theorem 1. By using Equation [28] and Equation |46j the joint density Ps^h^jV) °f tne 
ergodic invariant measure in terms of £ and i] is given as 

d(z, w, v) 



pz,H(£,,rj) 



pz,w,v(z,£,/z,v/ z ) 



dz 



Pz(z)p w ,v(£/z,ri/z)z 2 dz 



1 z exP2 -^(-r)2wwrw dz 

^ex P2 _ 1/? (--(C + r?) 



(50) 



Theorem 2. The marginal density of £ is a one-dimensional (/-Gaussian distribution, 



Pa(0 



' 1 - q' > 2 



1 f ^ 

exp ( — ^ 



i - ^Ke 

3-q'^ 



2tt 



i 1 i A (s-?') 



l + a^zif 2 



ki<a/^ ( 9 '<i) 



(q = q> = 1) 
^ (!<<?'< 3) 



(51) 



l-q' 2 ' 2 

where </ = (3q — l)/(q + 1). Hence, sequences generated from the maps in Definition 1. are random 
numbers sampled from a (/-Gaussian distribution, where q' = (3q — l)/(q + 1). 
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Proof of Theorem. 2 From Proposition 2., the marginal distribution of £ is the same functional form 
as Equation [T3l 



IV. Numerical simulation 



Figure [2] shows sample paths for several values of q'. As shown in these figures, they seem to be from 
a trapped random walk to Levy walk as q is increasing. Figure [3] shows the return maps between z n+ \ 
and z n . They show the determinism of the proposed random number generator. The return map of z n at 
I = 2 and c = 1 shows the functional form of the map function introduced in Equation [38] f2,i(z) = 
holds at z = J— 21n g (l/2). Since one has 



d/2,1 
dz 



2 1 -«(l-cxp g (-4))"'(l + (l-q)(-4)) T =? ; 
y / -21n,[2(l-exp,(-4))] 

(z< ^-21^(1/2)) 



2 1 -«2 



(52) 



-2]n,[2exp,(-£)] 

(z> ^-21^(1/2)) 



the Lyapunov exponent of z n , defined as 



1 



t-l 



n=0 



A = lim - V log f'(z n ) = log 2 = h KS , 

t— >oo t — 



(53) 



is computable. Here, hxs is the Kolmogrov-Sinai entropy. The relation A = Hks holds in one dimensional 
case by the Pesin identity. Independently of the initial conditions {vq,zq) and the parameter q, it is 
numerically confirmed that the Lyapunov exponent A approaches log(2) at I = 2 and c = 1. This is 
consistent with the theoretical value of chaotic map, which is conjugate with a diffeomorphism g for the 
tent map. More generally, the Lyapunov exponent A approaches to clog(Z) in a general case of fi ;C . This 
iterated map is deterministic, however, the auto-correlation function of the productive variable £ = wz, 



1 <_1 1 2 

C(m) = lim - €n£n+m - (lim - Y] £ n ) , 



(54) 



n=0 n=0 

decays for m > 1 from the orthogonality of the Chebyshev polynomials. Obviously, the expectation 
value of £ is 



1*1 ("OO 



(55) 
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Since due to the independence of w and z, we have 

1 

C(m) = lim - V w n z n P d o ■ ■ ■ o P d (w n ) fi >c o • • • o f l (z n ) 

n=0 Z. 

m in. 



w Pd o ■ ■ ■ o Pd(w) {w)dw 



z fi,c° • • • ° fi,c(z)p z {z)dz^J , 



(56) 



we obtain the auto-correlation of £ as 



C(m) = < 



l8 1 , 2m B{2, T ^) (q<l) 

^1,2- (9 = 1) 

|<5 1>2m B(2,a) (Kg<l) 



(57) 



Note that C(0) is not finite for 5/3 < q < 3 since the variance of g-Gaussian distribution is not finite 
for 5/3 < q < 2 and it is undefined for 2 < q < 3. In this derivation, we used the permutability and the 
orthogonality of the Chebyshev polynomials, 



/ Pa o • • • o Pj(w) P k o ■ ■ ■ o P k (w) fj,\y (w)dw 

J-\ > v ' V v ' 



Pj-n (w)P km {w)nw (w)dw 



~ jT COs(j n 0) COs(fc m fl)d0 = Itfj-n^ 



(58) 



In the same way, it can be proved that the auto-correlation function of the productive variable r/ = vz 
also decays for m > 1. 

The cumulative distribution of £ generated by Equation [37] Equation [38] and Equation [45] defined as 



Pr(~ <0 = [ Psit'W, 

J —CO 



can be expressed as 



Pr(~ < a = 



i erfc (A) 



3-<?' 

1-9' 



l + sign(C)/3(^' e • 1 1 



i I i'- 1 ei ' 2' o'-l 2/ 



for q = 1 
for 1 < q < 3 



(59) 



(60) 
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where f3(x;a,b) is the regularized incomplete beta function, 

P( x;a ,b) = — i— fV^l-t) 6 - 1 ^, (61) 
B{a,b) Jo 

and erfc(x) is the complementary error function denned as 

erfc(x) = — / e~* dt. (62) 

We compare the cumulative distributions of £ n obtained from Equation [37] Equation [38] and Equation 
|45l with Equation [60] Since we normally generate (/-Gaussian random variables from the given q', for 
practical usage, we need the inverse relation between q and q'\ q = (q' + l)/(3 — q'). Figure [4] shows the 
empirical complementary cumulative distributions of £, 

Pr(S>e) = l-Pr(3<0, (63) 

computed from 10,000 samples for (wq,zq) = (0.1, 1.0). Comparing the empirical distribution with the 
theoretical one, we found that they are very close for each parameter q'. 



We conducted the Kolmogorov-Smironov (KS) and the Anderson-Darling (AD) tests in order to verify 
whether the empirical distributions of sequences generated by our proposed method are convergent to the 
(/-Gaussian distributions. It is known that Anderson-Darling test is suitable for checking the goodness- 
of-fit for heavy-tailed distributions [22]. Assuming M samples of £i, . . . , £m> the test statistics are given 
as 

Z = vSmax 



F M (Cn) ~ Pr(H > £ n ) VV(Pr(£ > £«)), (64) 

where Fm{^u) an empirical cumulative distribution function, and t/j(u) is a weight function. In the case 
of ip(u) = 1, Z gives a KS test statistic and in the case of ip(u) = u ^i_ u ^ , Z gives an AD test statistic. 

Table U shows the best p-values of both KS and AD tests for several q values at d = 8, I = 2, and 
c = 1. The p-value of KS test is greater than 0.1 for q < 2.7. Therefore, the null hypothesis that the 
sequences are not samples from the theoretical distribution is not rejected at more than 5% statistical 
significance for q values from 1 to 2.6 in KS test. The degree of freedom v goes to as q approaches 3. 
For q > 2.7 (y < 0.17), both the proposed procedure and GBMM does not work since degree of freedom 
v is very small. The p- value of AD test is greater than 0.1 for q < 2.4. Since AD test is sensitive for tail 
events, the null hypothesis is not rejected from the value of q smaller than KS test values. Table ITTl shows 
the p-values of both KS and AD tests for several q values at d = 6, I = 2, and c = 6. The tendency 
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of p-values is very similar to ones at d = 8, I = 2, and c = 1. The KS test passes at more than 5% 
statistical significance for q values ranging from -1 to 2.6 in KS test. The same is true for — 1 < q < 2.4 
in the case of AD test. 

While GBMM fT"3l is based on transformation of uniform random variables, our proposed method here 
is purely mechanical generation of (/-Gaussian distribution based on ergodic theory. Thus, no random 
number are not assumed for the generations of (/-Gaussian distribution. Its implementation is very simple 
as shown in the example code in Appendix |Al Figures [6] (d = 8, / = 2, and c = 1) and [5] (d = 6, 
I = 2, and c = 6) show the best p-values of (a) KS test and (b) AD test obtained from 10,000 samples in 
100 trials with the proposal and the GBMM for several q. The best p-values provided by the proposed 
method are same as ones by the GBMM for many cases. 

V. Conclusion 

We proposed a pseudo random number generator of g-Gaussian random variables for a range of q values, 
— oo < q < 3, based on deterministic map dynamics. Our method consists of ergodic transformation on 
the unit circle and map dynamics based on the piecewise linear map. We conducted both KS and AD 
tests for random number sequences generated by GBMM and our proposed chaotic method for several 
values of q. The (/-Gaussian samples passed the KS test at the 5% significance level for q < 2.7, and 
passed the AD test at the 5% significance level for q < 2.4. 

Appendix A 
Source code 

We show a C source code for our proposed method for d = 8, I = 2, and c = 1. The code is exhibited 
in order to demonstrate the algorithm, and is not optimal for speed. The algorithm is implemented in four 
functions. The first two functions compute (/-exponential and (/-logarithmic functions. The next function 
setseed_qnormal(i;o, zq) sets two random seeds vq and zq, and qnormal(gO calls the iterated map to 
generate (/-Gauss random variables by our proposed method. 



#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
#include <strings.h> 

double qnormal_x, qnormal_y, qnormal_z; 
double expq (double q, double w) { 
if (q==1.0) { 

January 10, 2013 



DRAFT 



16 



return (exp (w) ) ; 

} 

else { 

return (expl (log (1 . 0+ (1 . 0-q) *w) / (1 . 0-q) ) ) ; 

} 

} 

double lnq (double q, double w) { 
if (q==1.0) { 

return (log ( w) ) ; 

} 

else { 

return ( (exp (log (w) * (1. 0-q) ) -1 . 0) / (1 . 0-q) ) ; 

} 

} 

void set seed_qnormal (double vO, double zO) { 
qnormal_x = sqrt ( 1-v0*v0 ) ; 
qnormal_y = vO; 
qnormal_z = zO; 

} 

double Q8 (double w, double v) { 

return (8*w*v* ( ( (16. 0*w*w-24 .0) *w*w+10 . 0) *w*w-l . 0) ) ; 

} 

double P8 (double w) { 

return ( ( ( (128. 0*w*w-256 . 0) *w*w+160 . 0) *w*w-32 . 0) *w*w+l . 0) ; 

} 

double f (double z){ 

return (1 . 0-fabs ( 1 . 0-1 . 99999*z) ) ; 

} 

void qnormal (double q) { 
double qq; 

qnormal_y = Q8 (qnormal_x, qnormal_y ) ; 
qnormal_x = P8 (qnormal_x) ; 
qq = (q+1 . 0) / (3 . 0-q) ; 

qnormal_z = f (expq (qq, -qnormal_z*qnormal_z*0 . 5 ) ) ; 
qnormal_z = sqrt (-2 . 0*lnq (qq, qnormal_z) ) ; 
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int main(int argc, char *argv[]){ 
double q, vO , z , eta, xi ; 
int i; 

if (argc != 4) { 

printf("%s q vO zO\n" , argv [ ] ) ; 
exit (0) ; 

} 

q = (double ) atof (argv [ 1 ]) ; 
vO = (double) atof (argv [2 ]) ; 
zO = (double) atof (argv [3] ) ; 
setseed_qnormal (vO, zO) ; 
for (i=0; K10000; i++) { 
qnormal (q) ; 

xi = qnormal_x*qnormal_z ; 
eta = qnormal_y*qnormal_z ; 
printf("%lf %lf\n", xi, eta) ; 

} 
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Fig. 2. A sample path of the map dynamics at d = 8, I = 2, and c = 1 for (a) q = -0.9, (b) -0.4, (c) 0.1, (d) 0.6, (e) 1.1 
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Fig. 3. Return map between z n and z„+i for (a) q' = -0.9, (b) -0.4, (c) 0.1, (d) 0.6, (e) 1.1 (y = 19), (f) 1.6 (v = 2.33), 

and (g) 2.1 (u — 0.818). The solid curve represents z„+i = h i(z n ) for each value of q — (q 1 + l)/(3 — q'). 
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Fig. 4. Complementary cumulative distribution functions of £ at d = 8, I = 2, and c = 1 for (a) q' = —0.9, (b) —0.4, (c) 0.1, 
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TABLE I 

The best KS and AD statistics obtained from 10,000 samples in 100 trials for several q at d = 8, 1 = 2, and 

C = 1. p- VALUES OF BOTH KS AND AD TESTS ARE SHOWN. 
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TABLE II 

The best KS and AD statistics obtained from 10,000 samples in 100 trials for several q at d = 6, 1 = 2, and 

C = 6. p- values of both KS AND AD tests are shown. 
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Fig. 5. (a) The best p-values of both (a) KS and (b) AD tests obtained from 10,000 samples in 100 trials with our proposed 
and GBMM for several q at d — 8, I — 2, and c = 1. 
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Fig. 6. (a) The best p-values of both (a) KS and (b) AD tests obtained from 10,000 samples in 100 trials with our proposed 
and GBMM for several q at d — 6, I — 2, and c = 6. 
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